

/*
 * Class:        MultinormalPCAGen
 * Description:  multivariate normal random variable generator
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */

package umontreal.iro.lecuyer.randvarmulti;

   import cern.colt.matrix.DoubleMatrix2D;
   import cern.colt.matrix.linalg.SingularValueDecomposition;

import umontreal.iro.lecuyer.probdist.NormalDist;
import umontreal.iro.lecuyer.randvar.NormalGen;
import umontreal.iro.lecuyer.rng.RandomStream;
import cern.colt.matrix.impl.DenseDoubleMatrix2D;


/**
 * Extends {@link MultinormalGen} for a <SPAN  CLASS="textit">multivariate normal</SPAN> distribution, generated via the method of principal components analysis
 * (PCA) of the covariance matrix. The covariance matrix 
 * <SPAN CLASS="MATH"><I><B>&Sigma;</B></I></SPAN> is 
 * decomposed (by the constructor) as 
 * <SPAN CLASS="MATH"><I><B>&Sigma;</B></I> = <B>V</B><I><B>&Lambda;</B></I><B>V</B><SUP>t</SUP></SPAN> where
 *  <SPAN CLASS="MATH"><B>V</B></SPAN> is an orthogonal matrix and 
 * <SPAN CLASS="MATH"><I><B>&Lambda;</B></I></SPAN> is the diagonal matrix made up
 * of the eigenvalues of 
 * <SPAN CLASS="MATH"><I><B>&Sigma;</B></I></SPAN>. <SPAN CLASS="MATH"><B>V</B><SUP>t</SUP></SPAN> is the transpose 
 * matrix of <SPAN CLASS="MATH"><B>V</B></SPAN>. The eigenvalues are ordered from the
 * largest (<SPAN CLASS="MATH"><I>&#955;</I><SUB>1</SUB></SPAN>) to the smallest (<SPAN CLASS="MATH"><I>&#955;</I><SUB>d</SUB></SPAN>). The random multinormal
 * vector 
 * <SPAN CLASS="MATH"><B>X</B></SPAN> is generated via
 * 
 * <P></P>
 * <DIV ALIGN="CENTER" CLASS="mathdisplay">
 * <B>X</B> = <I><B>&mu;</B></I> + <B>A</B><B>Z</B>,
 * </DIV><P></P>
 * where 
 * <SPAN CLASS="MATH"><B>A</B> = <B>V</B>()<SUP>1/2</SUP></SPAN>, and 
 * <SPAN CLASS="MATH"><B>Z</B></SPAN> is a <SPAN CLASS="MATH"><I>d</I></SPAN>-dimensional vector of
 * independent standard normal random variates. The decomposition method
 * uses the <TT>SingularValueDecomposition</TT> class in <TT>colt</TT>.
 * 
 */
public class MultinormalPCAGen extends MultinormalGen {
   private double[] lambda;

   private static SingularValueDecomposition getSvd (DoubleMatrix2D sigma) {
      return (new SingularValueDecomposition (sigma));
   }

   private DoubleMatrix2D decompPCA (SingularValueDecomposition svd) {
      DoubleMatrix2D D = svd.getS ();
      // Calculer la racine carree des valeurs propres
      for (int i = 0; i < D.rows(); ++i) {
         lambda[i] = D.getQuick (i, i);
         D.setQuick (i, i, Math.sqrt (D.getQuick (i, i)));
      }
      DoubleMatrix2D P = svd.getV();
      return P.zMult (D, null);
   }

   private void initL() {
      if (mu.length != sigma.rows() || mu.length != sigma.columns())
         throw new IllegalArgumentException
            ("Incompatible mean vector and covariance matrix");
      lambda = new double[mu.length];
      sqrtSigma = decompPCA (getSvd(sigma));
   }


   /**
    * Equivalent to
    *  {@link #MultinormalPCAGen((NormalGen, double[], DoubleMatrix2D)) MultinormalPCAGen}<TT>(gen1, mu, new DenseDoubleMatrix2D(sigma))</TT>.
    * 
    */
   public MultinormalPCAGen (NormalGen gen1, double[] mu, double[][] sigma) {
      super(gen1, mu, sigma);
      initL();
   }


   /**
    * Constructs a multinormal generator with mean vector <TT>mu</TT>
    *  and covariance matrix <TT>sigma</TT>. The mean vector must have the same
    *  length as the dimensions of the covariance matrix, which must be symmetric
    *  and positive semi-definite.
    *  If any of the above conditions is violated, an exception is thrown.
    *  The vector 
    * <SPAN CLASS="MATH"><B>Z</B></SPAN> is generated by calling <SPAN CLASS="MATH"><I>d</I></SPAN> times the generator <TT>gen1</TT>,
    *  which must be a <SPAN  CLASS="textit">standard normal</SPAN> 1-dimensional generator.
    * 
    * @param gen1 the one-dimensional generator
    * 
    *    @param mu the mean vector.
    * 
    *    @param sigma the covariance matrix.
    * 
    *    @exception NullPointerException if any argument is <TT>null</TT>.
    * 
    *    @exception IllegalArgumentException if the length of the mean
    *     vector is incompatible with the dimensions of the covariance matrix.
    * 
    */
   public MultinormalPCAGen (NormalGen gen1, double[] mu,
                             DoubleMatrix2D sigma) {
      super(gen1, mu, sigma);
      initL();
   }

   /**
    * Computes the decomposition <TT>sigma</TT> =
    * 
    * <SPAN CLASS="MATH"><I><B>&Sigma;</B></I> = <B>V</B><I><B>&Lambda;</B></I><B>V</B><SUP>t</SUP></SPAN>. Returns 
    * <SPAN CLASS="MATH"><B>A</B> = <B>V</B>()<SUP>1/2</SUP></SPAN>.
    * 
    */
   public static DoubleMatrix2D decompPCA (double[][] sigma)   {
      return decompPCA (new DenseDoubleMatrix2D (sigma));
   }


   /**
    * Computes the decomposition <TT>sigma</TT> = 
    * <SPAN CLASS="MATH"><I><B>&Sigma;</B></I> = <B>V</B><I><B>&Lambda;</B></I><B>V</B><SUP>t</SUP></SPAN>. Returns 
    * <SPAN CLASS="MATH"><B>A</B> = <B>V</B>()<SUP>1/2</SUP></SPAN>.
    * 
    */
   public static DoubleMatrix2D decompPCA (DoubleMatrix2D sigma)   {
      // L'objet SingularValueDecomposition permet de recuperer la matrice
      // des valeurs propres en ordre decroissant et celle des vecteurs propres de
      // sigma (pour une matrice symetrique et definie-positive seulement).

      SingularValueDecomposition sv = getSvd (sigma);
      // D contient les valeurs propres sur la diagonale
      DoubleMatrix2D D = sv.getS ();
      // Calculer la racine carree des valeurs propres
      for (int i = 0; i < D.rows(); ++i)
         D.setQuick (i, i, Math.sqrt (D.getQuick (i, i)));
      DoubleMatrix2D P = sv.getV();
      // Multiplier par la matrice orthogonale (ici P)
      return P.zMult (D, null);
   }


   /**
    * Returns the matrix 
    * <SPAN CLASS="MATH"><B>A</B> = <B>V</B>()<SUP>1/2</SUP></SPAN> of this object.
    * 
    * @return the PCA square root of the covariance matrix
    * 
    */
   public DoubleMatrix2D getPCADecompSigma() {
      return sqrtSigma.copy();
   }


   /**
    * Computes and returns the eigenvalues of <TT>sigma</TT> in
    *  decreasing order.
    * 
    */
   public static double[] getLambda (DoubleMatrix2D sigma) {
      SingularValueDecomposition sv = getSvd (sigma);
      DoubleMatrix2D D = sv.getS ();
      double[] lambd = new double[D.rows()];
      for (int i = 0; i < D.rows(); ++i)
         lambd[i] = D.getQuick (i, i);
      return lambd;
   }


   /**
    * Returns the eigenvalues of 
    * <SPAN CLASS="MATH"><I><B>&Sigma;</B></I></SPAN> in decreasing order.
    * 
    */
   public double[] getLambda() {
      return lambda;
   }


   /**
    * Sets the covariance matrix 
    * <SPAN CLASS="MATH"><I><B>&Sigma;</B></I></SPAN> of this multinormal generator
    *  to <TT>sigma</TT> (and recomputes 
    * <SPAN CLASS="MATH"><B>A</B></SPAN>).
    * 
    * @param sigma the new covariance matrix.
    * 
    *    @exception IllegalArgumentException if <TT>sigma</TT> has
    *     incorrect dimensions.
    * 
    * 
    */
   public void setSigma (DoubleMatrix2D sigma) {
      if (sigma.rows() != mu.length || sigma.columns() != mu.length)
         throw new IllegalArgumentException
            ("Invalid dimensions of covariance matrix");
      this.sigma.assign (sigma);
      this.sqrtSigma = decompPCA(getSvd (sigma));
   }


   /**
    * Generates a <SPAN CLASS="MATH"><I>d</I></SPAN>-dimensional vector from the multinormal 
    *  distribution with mean vector <TT>mu</TT> and covariance matrix
    *  <TT>sigma</TT>, using the one-dimensional normal generator <TT>gen1</TT> to
    *  generate the coordinates of 
    * <SPAN CLASS="MATH"><B>Z</B></SPAN>, and using the PCA decomposition of
    *  
    * <SPAN CLASS="MATH"><I><B>&Sigma;</B></I></SPAN>. The resulting vector is put into <TT>p</TT>.
    *  Note that this static method will be very slow for large dimensions, because
    *  it recomputes the singular value decomposition at every call. It is therefore
    *  recommended to use a <TT>MultinormalPCAGen</TT> object instead,
    *  if the method is to be called more than once.
    * 
    * @param p the array to be filled with the generated point.
    * 
    *    @exception IllegalArgumentException if the one-dimensional normal
    *     generator uses a normal distribution with <SPAN CLASS="MATH"><I>&#956;</I></SPAN> not equal to 0, or
    *     <SPAN CLASS="MATH"><I>&#963;</I></SPAN> not equal to 1.
    * 
    *    @exception IllegalArgumentException if the length of the mean
    *     vector is different from the dimensions of the covariance matrix,
    *     or if the covariance matrix is not symmetric and positive-definite.
    * 
    *    @exception NullPointerException if any argument is <TT>null</TT>.
    * 
    * 
    */
   public static void nextPoint (NormalGen gen1, double[] mu,
                                 DoubleMatrix2D sigma, double[] p) {
      if (gen1 == null)
         throw new NullPointerException ("gen1 is null");

      NormalDist dist = (NormalDist) gen1.getDistribution();
      if (dist.getMu() != 0.0)
         throw new IllegalArgumentException ("mu != 0");
      if (dist.getSigma() != 1.0)
         throw new IllegalArgumentException ("sigma != 1");

      if (mu.length != sigma.rows() ||
          mu.length != sigma.columns())
         throw new IllegalArgumentException
            ("Incompatible mean vector and covariance matrix dimensions");
      double[] temp = new double[mu.length];
      DoubleMatrix2D sqrtSigma = decompPCA(sigma);
      for (int i = 0; i < temp.length; i++) {
         temp[i] = gen1.nextDouble();
         if (temp[i] == Double.NEGATIVE_INFINITY)
            temp[i] = -MYINF;
         if (temp[i] == Double.POSITIVE_INFINITY)
            temp[i] = MYINF;
      }
      for (int i = 0; i < temp.length; i++) {
         p[i] = 0;
         for (int c = 0; c < temp.length; c++)
            p[i] += sqrtSigma.getQuick (i, c)*temp[c];
         p[i] += mu[i];
      }
   }


   /**
    * Equivalent to
    *  {@link #nextPoint((NormalGen, double[], DoubleMatrix2D, double[])) nextPoint}<TT>(gen1, mu, new DenseDoubleMatrix2D(sigma), p)</TT>.
    * 
    */
   public static void nextPoint (NormalGen gen1, double[] mu,
                                 double[][] sigma, double[] p) {
      nextPoint(gen1, mu, new DenseDoubleMatrix2D(sigma), p);
   }


   /**
    * Generates a point from this multinormal distribution. This is much
    * faster than the static method as it computes the singular value decomposition
    * matrix only once in the constructor.
    * 
    * @param p the array to be filled with the generated point
    * 
    */
   public void nextPoint (double[] p) {
      int n = mu.length;
      for (int i = 0; i < n; i++) {
         temp[i] = gen1.nextDouble();
         if (temp[i] == Double.NEGATIVE_INFINITY)
            temp[i] = -MYINF;
         if (temp[i] == Double.POSITIVE_INFINITY)
            temp[i] = MYINF;
      }
      for (int i = 0; i < n; i++) {
         p[i] = 0;
         for (int c = 0; c < n; c++)
            p[i] += sqrtSigma.getQuick (i, c)*temp[c];
         p[i] += mu[i];
      }
   }
}
